Expansion dynamics in the one-dimensional Fermi-Hubbard model 
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Expansion dynamics of interacting fermions in a lattice are simulated within the one-dimensional 
(ID) Hubbard model, using the essentially exact time-evolving block decimation (TEBD) method. In 
particular, the expansion of an initial band-insulator state is considered. We analyze the simulation 
results based on the dynamics of a two-site two-particle system, the so-called Hubbard dimer. Our 
findings describe essential features of a recent experiment on the expansion of a Fermi gas in a 
two-dimensional lattice. We show that the Hubbard-dimer dynamics, combined with a two-fluid 
model for the paired and non-paired components of the gas, gives an efficient description of the 
full dynamics. This should be useful for describing dynamical phenomena of strongly interacting 
Fermions in a lattice in general. 

PACS numbers: 71.10.Fd , 03.75.Ss, 73.20.Mf 



Important physical phenomena such as magnetism 
and high-temperature superconductivity are often ap- 
proached by theories based on the Hubbard model [H, Q 
which describes interacting particles in a lattice. Within 
ultracold gas systems [sl,!^? the Hubbard model can be ef- 
ficiently realized and studied in experiments with bosonic 
[El and recently with fermionic atoms [6, 7]. Intriguingly, 
the dimension can be easily controlled. Low-dimensional 
systems such as nanowires, iron pnictides and graphene 
are currently highlighted topics of research. Models for 
the quantum many-body physics of 2D and ID systems 
can explored with ultracold gases, c.f. recent experiments 
on fermions in one dimension [sl and expanding fermions 
in a 2D lattice [H*]. For one-dimensional systems, an ad- 
vantage is that the experiments can be compared to exact 
theoretical descriptions. However, although the ground 
state and static properties of one-dimensional systems are 
known to an impressive degree Q , dynamics is largely 
unexplored. Work on theory and simulation of dynami- 
cal properties of interacting fermions in ID has recently 
been emerging [12j . 

In this Letter, we study with exact numerical meth- 
ods the expansion of fermions within the one-dimensional 
Hubbard model. We show that the resulting complex dy- 
namics can be efficiently described by a two-fluid model 
in which we deduce the dynamics of the fluids from the 
dynamics of a Hubbard dimer. Our results explain sev- 
eral main features of the experiment [5] performed in 2D, 
and give exact predictions for future experiments in ID. 
The simple Hubbard-dimer two-fluid model that we have 
developed provides a basis for the description of various 
types of expansion, collision and oscillation dynamics for 
fermions in lattices. 

We use the time-evolving block decimation (TEBD) 
algorithm [l3| to describe the time evolution generated 
by the Hubbard Hamiltonian Hh = U hi^^hi^^ — 





FIG. 1: (color online) Schematic representation of the initial 
state: the middle part of the lattice is fully occupied {Oi) and 
the rest is empty {Ei). Sites El^ Ol and Or^ Er represent 
the left and right edge of the cloud, respectively. 



h.c.^ given an initial state |(/>(0)), 



where hi^a = Ci^Cia with and Cia representing the 
creation and annihilation of a fermion with spin a at 
the site i = 1 . . .L. Moreover, the initial state is given 

by \m) = I0)i • • • 10)^;. I t;)o. ... I t;)o«i0)s« . . . |0)l 

(see Fig. [1]). The initial state consists thus of a band in- 
sulator occupying the central Ol — Or sites of an other- 
wise empty lattice. In the simulation we have considered 
L = 150, El = 66, Ol = 67, Or = 86, Er = 87, the 
Schmidt number ^ = 150, and J = 1. Our code allows 
us to access the expectation value of the (spin-resolved) 
local particle number rii^it) and rtiiit) along with the lo- 
cal double occupancy Ui^iit) = {(j){t)\ni^ni i\(j){t)) . Note 
that nii{t) = ni^{t) since the problem is spin symmet- 
ric. In our analysis we will show that the evolution 
of the initial state can be described in terms of a two- 
fluid model where the two fluids are represented by sin- 
gle particles and douhlons as suggested by the Bethe- 
ansatz solution of the Hubbard model [3] and has been 
shown in the context of imbalanced Fermi gases 
18]. The doublons are excitations of the form cj^cj^|0) 

and the single (unpaired) particles are defined as 4^|0) 
(cr =t, I). The local number of doublons is given by 
ni^l(t), while the number of unpaired (up) particles is 
given by n^^(t) = ni^{t) - ni^i{t). 

In Figs. [2] and [3l ^/ni^{t) and are depicted 

for = 0.0 and U = ±10.0. We are plotting the square 
roots of the densities since this highlights low density 
features (see the supplementary material for the full den- 
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FIG. 2: (color online) Time evolution of y^ni^{t) (above) and 
y^ni^l{t) (below) for \U\ = 0.0. The free-particle nature of 
the expansion is clear from the absence of separate doublon 
expansion wavefronts. 
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FIG. 3: (color online) Same as Fig. [2]for \U\ = 10.0. For this 
interaction it is possible to distinguish the two wavefronts. 



sity plots). As in general for the spin-balanced Hubbard 
model [3] , the density distributions evolve in time exactly 
in the same way for U and —U. This U ^ —U symme- 
try holds also for all observables for all interactions and 
it was indeed observed also in the experiment [5]. 

In the non-interacting case in Fig. [21 both particles and 
doublons are expanding at the speed of 2 J, corresponding 
to the highest group velocity allowed by the dispersion 
relation (see supplementary material). For strong inter- 
actions (Fig. [3]), we see two separate wavefronts. Such 
separation into two types of wavefronts is clearly observ- 
able for interactions \U\ > 3.0. The outermost wavefront 
consists of fully unpaired particles expanding at the of 
speed 2J, like in the U = case. In contrast doublons 
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FIG. 4: (color online) Unpaired particle expansion n^^(t) for 
\U\ = 5.0, exhibiting the oscillations described in the text. 
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FIG. 5: Unpaired population dynamics n^^^^^ ^(t) for 
\U\ = 5.0 and \U\ = 10.0. 



expand at the speed of AJ^ /U (see supplementary mate- 
rial for the explanation of these results). 

Intermediate interactions 0.5 < \U\ < 3.0 show a more 
complicated behaviour. The separate expansion fronts 
are no longer well distinguishable, suggesting a stronger 
interplay between single particles and doublons. Mov- 
ing to even lower interactions, the unpaired particles and 
the doublons behave similarly to the case of U = 0. 
Let us now examine the time dependence of the num- 
ber of unpaired (up) particles n^^{t) for \U\ = 5.0 (Fig. 
Sj). Initially the dynamics in the band insulator cloud is 
Pauli blocked, since neighbouring lattice sites in the cen- 
ter have unit density for both spin up and down. There- 
fore the unpaired expansion fronts are created at the 
edges of the cloud. Intriguingly, n^^it) shows damped 
oscillations at the edges associated with the emission of 
unpaired particles into the empty lattice. Considering 
the time evolution of ri^^ ^ + n^J: ^ (the two edge sites) 



over the whole interaction range \U\ = 0.0 — 15.0 (see 
Fig. [5] for \U\ = 5.0, 10) we find that there are fewer pe- 
riods of oscillations for lower interactions (for \U\ = 1.0 
only one broad oscillation peak is visible) and the oscil- 
lation frequency increases with interaction strength, see 
supplementary material for a general survey of the data. 

The evolution of doublons into unpaired particles plays 
a key role in the expansion physics. For this reason, we 
focus now on the explanation of the oscillations in the 
case of high interactions, see Figs. [4land[5l Our hypoth- 
esis is that one can consider the edges of the cloud (sites 
El^ Ol and Or^ Er) at short timescales as two-site sys- 
tems (Hubbard dimers [1, Focusing on the El/Ol 
dimer, the system can be described as an initially empty 
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FIG. 6: The frequency given by the FT of n^'^+Oi.,t(^) (^o^^ 
TEBD numerics) compared to the frequency ^JV^ + 16J^ ob- 
tained by solving the two-site model. 
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FIG. 7: The height of the first peak of nl;^+o,.,t(^) (^^^ 
TEBD numerics) compared to the amplitude obtained by 
solving the two-site model. 

state |0) in the left lattice site Ei^ and a doublon | in 
the right lattice site Ol- The dynamics of the dimer with 
this initial state can be solved analytically by diagonal- 
izing its Hamiltonian (see supplementary material). As 
a result, in the two-site problem, the population of un- 
paired up particles on the two sites Ej, and Ol is given 
by (the tilde refers to the Hubbard dimer model): 

_ l-cos(V?7^ + 16J^t) 

^ + 8J2 

We extract the oscillation frequencies from the numer- 
ical Fourier transform (FT) of ^l;^+Oi,,t compare 
its peaks to the oscillation frequency \JU^ -h 16J^ in ([T]) 
(see Fig. [6]). Moreover, we compare the height of the first 
peak of the unpaired density oscillations (seen in Figs. S] 

andO to the amplitude 8/ (^16 + in Fig. [71 

The agreement is good considering that, for longer 
times, the FT of ^^^+o^ ^(^) has additional contribu- 
tions stemming from the hopping between the dimer and 
the rest of the chain. However, as the frequency is ap- 
proximately U in the high interaction limit, one might 
claim that the frequency correspondence could be ob- 
tained from simple energy arguments. Therefore, to fur- 
ther confirm the validity of the Hubbard Dimer model, 
we now move on to examine whether the two-site model 
coupled to the next adjacent sites explains the decay ob- 
served in ^l;^+Oi,,t(^) ^^^^ define the 
damping D to mean the decrease of the amplitude of the 
unpaired density oscillations, compared to t = 0. We pro- 



pose that the damping should be equal to the probabil- 
ity of having an unpaired particle in the Hubbard dimer 
times the probability for this single particle tunnelling 
out of this system. The probability for the single particle 
tunnelling (obtained by solving the two-site system with 
the initial state |0,t)) is given by sin^(Jt). Combining 
this result with Eq.([T]), we obtain the damping at a given 
time r: 

N r 1 - COs(V/72 + \UH) . 2/ T N . 

D(t) = 2 / "-^-jj, ^ sm\jt)dt, (2) 

where the factor of two comes from the particle-hole sym- 
metry: particles leaking out from the dimer to the left 
(to El — 1) are mirrored by holes leaking to the right 
(to Ol + 1) thus generating particle-like expansion fronts 
emitted out of the initial cloud and hole-like expansion 
fronts emitted into the cloud. When the particle-like and 
hole-like expansion fronts meet interference in the un- 
paired particle density is visible, see Fig. |4j 

By comparing the decay predicted in Eq. (|2]) to the 
numerics, we observe that, for the duration of the first 
half-period, the decay is negligible. This is in accordance 
with the height of the first peak being equal to the two- 
site oscillation amplitude, as shown in Fig. [71 After three 
half-periods we compare the damping as predicted by Eq. 
([2|) to the change of amplitude between the first peak and 
the second peak seen in numerics, see Fig. [H The two- 
site model is again in good agreement with TEBD for 
\U\ > 3.0. The time beyond which Eq. (j2j) fails to de- 
scribe the expansion physics is when the population of 
unpaired particles in the sites Ol+i,^l-i becomes sig- 
nificant. This occurs when the siii?{Jt) = ^~^^^^^(^-^^) 
term is no longer close to zero, limiting our short time 
analysis to times t « ^ j. When U is sufficiently large, 
the Hubbard dimer oscillations occur in much shorter 
timescale than the single particle tunneling does. In other 
words, a large number of oscillations occur in the win- 
dow < t < ^ J In the case of lower interactions, Hub- 
bard Dimer oscillations at frequency VW^ + 16J^ become 
comparable to the frequency 2 J and therefore we see that 
already the first oscillation peaks are heavily damped. 

In general, the two-site dynamics is well able to de- 
scribe the creation of the particle, hole, and doublon 
wave fronts seen in the density profiles. These wave fronts 
are created during the two-site oscillations. Our numeri- 
cal results and analytical investigations confirm the two- 
fluid picture of the system. Initially, the interaction takes 
place at the edges of the cloud, where unpaired parti- 
cles are created according to the dimer dynamics de- 
scribed above. The subsequent expansion is explained 
by dynamics of non-interacting particles (at the speed of 
2 J) or doublons (at the speed of 4J^/t/). Our model 
gives an excellent quantitative description in the highly- 
interacting limit due to the clear separation of the ex- 
pansion and dimer oscillation eigenfrequencies. For in- 
teractions 0.5 < \U\ < 3.0 the interplay between the 
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FIG. 8: Change of the amphtude of n^^^^(t) oscillations after 
t — ' compared to the change of amplitude between 

the first and second oscillation peaks of ^o^+£;^,^(t) observed 
in TEBD numerics. 
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FIG. 9: The core expansion speed i?(t), as a function of the 
interaction strength. For reference, we plot also the expansion 
speed of the whole cloud. 



expansion and the Hubbard dimer dynamics does not al- 
low a quantitative description of the numerical results, 
however it provides a qualitative framework for further 
analysis. For \U\ < 0.5, the free-particle expansion seems 
to give a fairly good description. 

Finally, we compare our results to the 2D experiment 
of [5] done at finite temperature. We suggest that our 
two-site considerations also apply to the dynamics of the 
experiment. In 2D there is coupling to four adjacent 
sites, but just like in ID, initially the sites are either Pauli 
blocked or empty. The simplified dynamics should orig- 
inate from the two-site analysis, and subsequent short- 
time dynamics in the high interaction limit correspond 
to the two-fluid expansion, with the two fluids interact- 
ing only at the edges. 

In order to compare our results to Fig. 5 of [5] we 
define the core density as np^(t) = ni^(t) — n\^{t) for 
\U\ > 0.5 and nf^{t) = ni^{t) for \U\ < 0.5, where n^^^{t) 
is the density of purely ballistic particles (see supplemen- 
tary material). The definition of the core then corre- 
sponds to the diffusive part of the cloud in the model 
of [Bl]. The core expansion velocity is given by R{t)^ 
where R{t) = v'<72^^^^T>2, < i >= Ezti(k^tW + 

^u(^)] * ^) / Ylf=ii'^?ti'^) + ^u(^)) ^ number 
of lattice sites. The core expansion speed as a function 
of interaction is plotted in Fig. [9l The behavior of R{t) 
is indeed similar to the core expansion velocity in Fig. 5 
of [H, showing also the negative velocities. 

Another recent experiment [loj studied collision dy- 
namics of two Fermi gas clouds. Although the experi- 
ment is not done in a lattice, the theoretical framework 
presented here can in the low density limit be used to 
desribe also physics in [lo|, c.f. [loj . 

In conclusion, we studied the expansion of an interact- 
ing fermionic gas in a ID lattice. We showed that the 
time evolution of this system can be described in terms 
of a two-fiuid model of unpaired particles and doublons 
whose interplay gives rise to nontrivial dynamics. We 
suggest that the experimental results of [5] can be in- 
terpreted in terms of the analysis performed here. Our 
results should be widely applicable since expansion is a 



basic dynamics problem related to the ultracold gas ex- 
periments in particular, and to the transport properties 
of fermions in general. 

We thank I. Bloch and U. Schneider for useful dis- 
cussions. This work was supported by the Academy of 
Finland (Projects No. 213362, No. 217043, No. 217045, 
No. 210953, and No. 135000) and EuroQUAM/FerMix, 
and conducted as a part of a EURYI scheme grant (see 
www.esf.org/euryi). The research was partly supported 
by the National Science Foundation under Grant No. 
PHY05-51164. Computing resources were provided by 
CSC the Finnish IT Centre for Science. 



* Electronic address: paivi.torma@hut.fi 
[3] F. H. L. Essler, H. Frahm, F. Gohmann, A. Kliimper, 

and V. E. Korepin, The One- Dimensional Hubbard Model 

(Cambridge University Press, 2005). 
[2] K. L. Hur and T. M. Rice, Ann Phys-New York 324, 1452 

(2009). 

[3] D. Jaksch and P. Zoller, Ann Phys-New York 315, 52 
(2005). 

[4] I. Bloch, Nature Physics 1, 23 (2005). 

[5] M. Greiner, et al.. Nature 415, 39 (2002). 

[6] R. Joerdens, et al.. Nature 455, 204 (2008). 

[7] U. Schneider, et al.. Science 322, 1520 (2008). 

[8] Y. an-Liao, et al.. Nature 467, 567 (2010). 

[5] U. Schneider, et al., arXiv cond-mat. quant-gas 
1005.3545vl, (2010). 

[10] A. Sommer, et al., arXiv cond-mat. quant-gas 
1101.0780vl, (2010). 

[2] T. Giamarchi, Quantum Physics in One Dimension (Ox- 
ford University Press, 2003). 

[12] C. Kollath, U. Schollwock, and W. Zwerger, Phys. Rev. 
Lett. 95, 176401 (2005); F. Heidrich-Meisner, et al., Phys. 
Rev. A 78, 013620 (2008); F. Massel, M. J. Leskinen, 
and P. Torma, Phys. Rev. Lett. 103, 066404 (2009); 
F. Heidrich-Meisner, et al., Phys. Rev. A 80, 041603 
(2009); A. Yamamoto, M. Yamashita, and N. Kawakami, J 
Phys Soc Jpn 78, 123002 (2009); M. Tezuka and M. Ueda, 
New J Phys 12, 055029 (2010); A. Korolyuk, F. Mas- 
sel, and P. Torma, Phys. Rev. Lett. 104, 236402 (2010); 



5 



S. Diehl, et aL, Phys. Rev. Lett. 105, 227001 (2010). 
[13] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); A. Daley, 

et al., J Stat Mech-Theory E p. P04005 (2004). 
[14] E. Zhao and W. V. Liu, J Low Temp Phys 158, 36 (2010). 
[15] G. Orso, Phys. Rev. Lett. 98, 070402 (2007). 
[16] H. Hu, X.-J. Liu, and P. D. Drummond, Phys. Rev. Lett. 

98, 070403 (2007). 
[17] X.J. Liu, H. Hu and P.D. Drummond, Phys. Rev. A 76, 

043605 (2007). 

[18] X.W. Guan, M.T. Batchelor, C. Lee, and M. Bortz, Phys. 

Rev. B 76, 085120 (2007). 
[19] J. Kajala, F. Massel, and P. Torma, arXiv cond- 

mat.quant-gas 1102.0133, (2011). 



EXPANSION DYNAMICS IN THE 
ONE-DIMENSIONAL FERMI-HUBBARD 
MODEL - SUPPLEMENTARY MATERIAL 

SURVEY OF THE NUMERICAL RESULTS 

Density profiles 

We show in Figs. [TOl and [TTI the time evolution of the 
(square-root) density for unpaired particles ^Jni^{t) and 
doublons ^Jni^l{t) for different values of the interaction. 
While the outermost wavefront velocity (2J), being re- 
lated to the expansion of the unpaired particles, in Fig. 
[To] is interaction-independent, the doublon wavefronts in 
Fig. [TT]have an interaction-dependent velocity (4J^//7). 
In order to give grounds for our choice to plot the square 
root of unpaired particles and doublons in the article, 
in Figs. [12] and [13] we compare the square root densi- 
ties and the densities for \U\ = 10.0. It is noted that 
we have studied the time evolution up to time t = 30 
1/J. At this time the expanding cloud has not reached 
the edges of our finite system. However, the finite size 
causes low- momentum cutoff in the expansion dynamics. 
Therefore, we have experimented with systems of differ- 
ent size and we have seen that the lattice size does not 
affect the expansion dynamics significantly 



Doublon and unpaired particle oscillations 

In the main article, we compare the oscillation of the 
unpaired population in the Hubbard dimer ^l;^+Oi,,t(^) 
to its counterpart obtained with the TEBD algorithm 
^l;^+OL,t(^)* [13]^^^ number of doublons in the 

whole lattice nj^^]^\t) as a function of time for various 
interactions are shown, exhibiting increasing oscillation 
frequency with increasing \U\. In Figs. [T5] and [T6] we 
show the Fourier transform (FT) of TiEL+OL,^i{^) 
n^ll^^(t) for \U\ = 5. Owing to particle- number con- 
servation ^l;^+Oi,,t(^) ^^o^^ same time oscillatory 
dependence as TiEL+OL,^i{^)i although naturally in the 
opposite phase. The presence of extra peaks present 



in the FT of TiEL+OL,^i{^) with respect to the FT of 
^^l+Ol t^^) ^^^^ explained later. 
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FIG. 10: The square root of up density profiles y^ni^{t). 
\U\ = 0, 0.1, 1.0, 5.0, 10.0 (top left to bottom right) 
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FIG. 11: The square root of doublon density profiles. \U\ 
0, 0.1, 1.0, 5.0, 10.0 (left to right, top to bottom) 
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FIG. 12: Comp arison between the square root of the density 
profile y^ni^{t) (left) with the density profile ni^(t) (right) for 
|[/| = 10.0. The unpaired particle wavefronts expanding with 
speed 2 J are clearly visible in the former. 
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FIG. 13: Comparison between the square root of the dou- 
blons density profile ^/nTuJt) (left) with the density profile 
^*n(t) (right) for \U\ = 10.0. In the left panel ( y^m^it)) 
low-density features are enhanced. 
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FIG. 14: The time dependence of n^tf^(^) 1^1 = 
0.1, 1.0, 5.0, 10.0 (top left to bottom right) 




-15 -10 -5 5 10 15 

Angular frequency (J) 



FIG. 15: The Fourier transform of n£;^+o^,ti(^) for U=5.0. 



In the numerics the FT n'^^_^Q^ ti^^) shows extra 
peaks with respect to the dimer dynamics which are asso- 
ciated to the hopping between the dimer and the rest of 
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FIG. 16: Fourier transform of n^f ^^(t) for U=3.0 (left) and 
U=5.0 (right) 
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FIG. 17: The position of the peaks in the Fourier transform 
of riE]^ as a function of the interaction strength. 



the chain. The peak at :^ [/, present both for n^^^^\t) 
and TiEL+OLAi{^)^ correspond to the doublon dissocia- 
tion frequency, as it will be later shown. 



THE HUBBARD DIMER MODEL 

The eigenstates of the two-site Hubbard Hamiltonian 
(so called Hubbard dimer) 

H = Hj + Hint 

Hj = -J^^c\^^C2a + h.C. 
a 

Hint = U ^ Ui^riii (3) 

i=l,2 
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are given by (see e.g. 



|T) = — (|t,i) + U,t)) 

\D.) = ^{\n,o)-\o,n)) 



\v±) = 



1 



l + al 



(|5)+a±|D+)), 



(4) 



where 



a± 



U_ 
47 



1±\ 1 



16J2 

C/2 



l5) = -^(it,;)-u,t)) 
l^+> = ^(in,o) + |o,n» 



(5) 
(6) 
(7) 



and the corresponding eigenvalues have the following 
form 



A_ = U/2 

Ao =0 
Xu = U 

A+ = U/2 



16J2 



16J2 



(< 0) 



(> u) 



\v-) 

\T) 

\D-) 



(8) 



Hubbard dimer dynamics 



We can apply the approach described in the previous 
section to the analysis of the edge sites of the cloud. As 
described in the main article, the edges of the cloud can 
be described as Hubbard dimers. If we focus on the sys- 
tem El -\- Ol^ the initial state is given by 

m = 0)) = |0,n> = ;^(7+l«+) - l-\v-) - \D-)) (9) 

where we have denoted 



7± 



(10) 



The number of doublons in the initially empty site El is 
given by 



(0,n|e"^n^;e-V|0,n) 

4 ^ 4(q!+-q;_)2 ^ 2{a+-a-y 
q;_)_q;_ +1 



COS [(A+ - X-)t] 
cos [{Xu - A+)t] 



+ 2(a;:a-) COs[(Aa-A_)t]. (11) 



Straightforwardly, Eq. (pTj) allows to evaluate the number 
of unpaired particles in the same site 



OtJ^CX- 



2(a+-a_)2 
1 



"2( 



2{a 



2(a+-Q!_)2 ^ 2(q!+ 
*C05 [(A+ — \-)t] 

y C0S[(A[/ - A+)t] 

COS [(Af/ — A_)t] . 



-) 



(12) 



The extra peaks appearing in the FT of riEL+OL^^iif) 
(Fig. [IS]), with respect to the FT of n^f ^^(t) (Fig. [IS] 
are due to the fact that, on top of the doublon dissoci- 
ation dynamics, ^£;i,+OL,ti(^) affected by the hopping 
of both unpaired particles and doublons in and out of the 
dimer El/Ol- These contributions, since a sum over the 
whole cloud is carried out, add up to zero for ri^^^^\t). 



EXPANSION VELOCITIES FOR PARTICLES 
AND DOUBLONS 

The expansion velocity of the unpaired particles can 
be easily obtained considering the dispersion relation of 
a particle hopping through on an empty lattice 



e/c = —2Jcos{k) 



leading to the group velocity 



^/c ^ = 2J sin(/c). 
ok 



(13) 



(14) 



Since in our system we consider particles generated at a 
specific site (El)^ all momentum states contribute to the 
expression of the initial state. From Eq. ([T4]) we thus see 
that the maximum velocity is attained for /c = 7r/2 



2J 



(15) 



which corresponds to the value obtained in the TEBD 
simulation. Analogously it can be shown that, for our 
initial configuration, the doublon hopping corresponds 
to the hopping of a particle with Jeff = ^ [l-Q cor- 
responding to the velocity of propagation of a charge ex- 
citation in the strongly attractive Hubbard model. In 
particular, if we consider the case where no unpaired 
particles are present, we can map our system to an 
isotropic Heisenberg chain by identifying | ti) =^ I t) 
and |0) =^ I ^). To prove this, we consider: 



1. A staggered particle- hole transformation 



(t) 



(-i)Vu 

(-1)^4 



i 4- 



(16) 
(17) 
(18) 
(19) 
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which transforms our original system (fihing n = 
^.^n^cr/^, magnetization density m = and in- 
teraction strength [/) to a half-fihed system with 
fining = 1, magnetization density = n — 1 and 
interaction strength U' = —U . Since n'^^ = ^i^i 

2. by strong coupling second order expansion the 
system Hamiltonian is then mapped to an isotropic 
Heisenberg Hamiltonian 



Hh 



u 



(20) 



With the mappings described above, since a doubly 
occupied site correspond to a spin up, and an empty 
site to a spin down, the band-insulator initial state 
10(0)) is thus transformed to a spin chain with Or — 
Ol central spin- up sites surrounded by spin-down 
sites, and the expansion dynamics of the doublons 
corresponds to the propagation of a single spin-up 
excitation through a spin-down polarized system. 
Its dispersion relation is given by 



— cos{k) 



(21) 



and hence, analogously to the derivation of the 
maximum velocity for the unpaired particles, to 



V 



doublon ^spin 



MAX 



MAX 



4J2 



U 



(22) 



DEFINITION OF THE CORE 

We have in the high interaction limit ballistic particles 
emitted from the diffusive core, the two types of ballistic 
particles being unpaired particles and doublons. These 
two expand at constant speeds (2 J and ^) in wave- 
fronts whose densities are so low that spin-spin scatter- 
ing (i.e. the Hubbard Dimer dynamics) is negligible. In 
contrast, the core has high density, the Dimer dynamics 
occur, and thus the core is diffusive. In the case of in- 
teractions lower than \U\ < 3 also the densities of the 
unpaired and doublon wavefronts become high enough 
for the Dimer Dynamics, which in the language of the 
classical non-linear diffusion model in [5] makes the two 
types of wavefronts diffusive. This diffusive behaviour 
is seen in our results as the absence of ballistic wave- 
fronts in density profiles (see e.g. the \U\ = 1.0 case 
in Fig. [To] above). Finally, for the very low interactions 
U < 0.5 the problem simplifies to ballistic expansion of 
non- inter acting particles. 

In order to compare our results to Fig. 5 of [5] we 
arrive at the problem of defining the core. In the high 
interaction limit this is perhaps easier as the ballistic par- 
ticles clearly separate from the central diffusive region. 





- Expansion speed after removing 

all ballistic contributions 
-A- Expansion speed after removing 

the ballistic particle-singlets 
■ ■ + ■ Expansion speed of the whole cloud 











5 10 
Interaction |U| (J) 



15 



FIG. 18: The core expansion velocities obtained by removing 
ballistic contributions (see text). 



the core. In the middle interaction region (see again the 
U = 1.0 case in Figs. [10] and [11] below ) we see that 
only one outermost wavefront expands ballistically and 
therefore if we consider the diffusive region of the cloud 
to be the core then the core would be everything else 
but the outermost wavefront. This definition however is 
problematic in the low interaction limit (U < 0.5). In 
this limit everything can be described as expansion of 
non-interacting ballistic particles. So the definition of 
the core as the diffusive part of the cloud does not work 
any more as then we would not have a core at all. 

However, experimentally it might be easier to define 
a high density region to correspond to the core, and fit 
the half-power half-maximum of a Gaussian to give the 
core width and subsequently the core expansion speed. 
In effect, then, in the low interaction regime the core is 
defined as the whole cloud, and in the middle and high 
interaction regimes the core is defined as everything else 
except the ballistic expansion fronts. This is what we 
have done in the main article. 

A variation of the above way of removing ballistic con- 
tributions would be to remove only the ballistic particle 
wavefronts, as the ballistic hole wavefronts may in effect 
be included in the experimental high density core region. 
To elaborate on what the hole wavefronts are, in the core 
state I ti) changes into a superposition of | t> and | 1). 
This decreases the total density, but in addition the | t) 
and I D states created are now part of a singlet state 
which expands ballistically into the core. So, if we want 
to define the core to be everything that is not ballistic, we 
need also to remove the density contribution of the bal- 
listic singlets in the core. Thus there are two reasonable 
approaches in defining the core. First is that we remove 
everything that is ballistic (particles and holes) and the 
second is that we remove only the particles. The core 
expansion speeds obtained in these two ways are plotted 
in [18] Both methods result in negative core expansion 
velocities as seen in the experiment [5i]. 
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THE TIME EVOLUTION OF CHARGE-CHARGE 
CORRELATIONS 

Using TEBD we can calculate correlation functions 
during the expansion of the cloud. An interesting correla- 
tion function in the case of the band insulator expansion 
is the charge-charge correlator given by 

a,j{t) =< $(t)|(n,,t +na)Kt + > • (23) 

However, Cij{t) is dominated by single particle den- 
sity correlations. To see the effect of the non- trivial cor- 
relations, we want to remove the density-density corre- 
lations from Cij{t). If i ^ j the density-density corre- 
lations are given by < hi^^ 

> + < ^i,t >< 
When i = j, < rii^^rii^^ > 



>< 
> 4 



- < Ui^l >< 



Ci,t > = < c. 



(1 



>=< cl^Ci^^ >= rii^^ and therefore on the 
diagonal the we just remove the single particle and dou- 
blon densities. Therefore, the correlation function which 
displays the non-trivial correlations is 



(< hi^^ >< hj^^ > 
+ < ni,; >< nj,; > 
< hi^^ >< fij^^ > 
+ < ni,; >< fij^^ >) 

-(< ni,^ > + < hi^i > 
+2 < fii^^i 



(24) 



The square root of the correlator Dij{t) is plotted in 
Figs, [inland [20] for \U\ = 5.0 at times 4^,8^,12^,16^. 
At time t = Oj the correlator is zero everywhere. We plot 
the square root instead of the actual value to see better 
the low value regions of the correlator. While keeping 
r = 150 during the time evolution in order to minimize 
the truncation error propagation, for the calculation of 
the correlators we have truncated the Schmidt number 
to r = 80 due to heavy numerical cost of the calculation, 
since for the evaluation of an observable at a given time 
the introduced error is negligible. Our dimer picture is 
supported by the fact that most of the correlations arise 
between sites around the edges of the cloud. In addition 
to this, correlations arise between the wavefronts expand- 
ing in the initially free sites i < El and i > Er and the 
whole initial cloud, corresponding to the stripes appear- 
ing in Figs. [T9land[2Ql A similar description can be ap- 
plied to the holes propagating towards the center of the 
cloud. This behavior reflects the intrinsically collective 
nature of the excitations in ID systems. Our picture of a 
ballistic particle leaving the cloud and a hole propagating 
towards the centre is complemented by this observation, 
which leads to the buildup of correlations between the 
propagating particles and the entire cloud. 





FIG. 19: (color online) Left: y^D~(t) for \U\ = 5.0 at t = 
4^. Right: The same for t = S^^. 
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FIG. 20: (color online) Left: The same as figure [19] but for 



t=12j (left) and t : 



16j 



jight). 
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